function [Pxx, Pxxc, f] = pspecd(varargin)
%PSD Power Spectral Density estimate.
%   Pxx = PSD(X,NFFT,Fs,WINDOW) estimates the Power Spectral Density of 
%   signal vector X using Welch's averaged periodogram method.  X is 
%   divided into overlapping sections, each of which is detrended, then
%   windowed by the WINDOW parameter, then zero-padded to length NFFT.  
%   The magnitude squared of the length NFFT DFTs of the sections are 
%   averaged to form Pxx.  Pxx is length NFFT/2+1 for NFFT even, (NFFT+1)/2
%   for NFFT odd, or NFFT if the signal X is complex.  If you specify a 
%   scalar for WINDOW, a Hanning window of that length is used.  Fs is the 
%   sampling frequency which doesn't affect the spectrum estimate but is 
%   used for scaling of plots.
%
%   [Pxx,F] = PSD(X,NFFT,Fs,WINDOW,NOVERLAP) returns a vector of frequen-
%   cies the same size as Pxx at which the PSD is estimated, and overlaps
%   the sections of X by NOVERLAP samples.
%
%   [Pxx, Pxxc, F] = PSD(X,NFFT,Fs,WINDOW,NOVERLAP,P) where P is a scalar
%   between 0 and 1, returns the P*100% confidence interval for Pxx.
%
%   PSD(X,...,DFLAG), where DFLAG can be 'linear', 'mean' or 'none', 
%   specifies a detrending mode for the prewindowed sections of X.
%   DFLAG can take the place of any parameter in the parameter list
%   (besides X) as long as it is last, e.g. PSD(X,'mean');
%   
%   PSD with no output arguments plots the PSD in the current figure window,
%   with confidence intervals if you provide the P parameter.
%
%   The default values for the parameters are NFFT = 256 (or LENGTH(X),
%   whichever is smaller), NOVERLAP = 0, WINDOW = HANNING(NFFT), Fs = 2, 
%   P = .95, and DFLAG = 'none'.  You can obtain a default parameter by 
%   leaving it off or inserting an empty matrix [], e.g. PSD(X,[],10000).
%
%   See also CSD, COHERE, TFE, PMEM, PMTM, PMUSIC.
%   ETFE, SPA, and ARX in the Identification Toolbox.

%   Author(s): T. Krauss, 3-26-93
%   Copyright (c) 1988-98 by The MathWorks, Inc.
%   $Revision: 1.27 $  $Date: 1997/12/02 18:37:07 $

%   The units on the power spectra Pxx and Pyy are such that, using
%   Parseval's theorem: 
%
%        SUM(Pxx)/LENGTH(Pxx) = SUM(X.^2)/LENGTH(X) = COV(X)
%
%   The RMS value of the signal is the square root of this.
%   If the input signal is in Volts as a function of time, then
%   the units on Pxx are Volts^2*seconds = Volt^2/Hz.
%
%   Here are the covariance, RMS, and spectral amplitude values of
%   some common functions:
%         Function   Cov=SUM(Pxx)/LENGTH(Pxx)   RMS        Pxx
%         a*sin(w*t)        a^2/2            a/sqrt(2)   a^2*LENGTH(Pxx)/4
%Normal:  a*rand(t)         a^2              a           a^2
%Uniform: a*rand(t)         a^2/12           a/sqrt(12)  a^2/12
%   
%   For example, a pure sine wave with amplitude A has an RMS value
%   of A/sqrt(2), so A = SQRT(2*SUM(Pxx)/LENGTH(Pxx)).
%
%   See Page 556, A.V. Oppenheim and R.W. Schafer, Digital Signal
%   Processing, Prentice-Hall, 1975.

error(nargchk(1,7,nargin))
x = varargin{1};
[msg,nfft,Fs,window,noverlap,p,dflag]=psdchk(varargin(2:end),x);
error(msg)

% compute PSD
window = window(:);
n = length(x);		% Number of data points
nwind = length(window); % length of window
if n < nwind    % zero-pad x if it has length less than the window length
    x(nwind)=0;  n=nwind;
end
x = x(:);		% Make sure x is a column vector
                        % this line must be AFTER the zero-padding in
                        % case x is a scalar
k = fix((n-noverlap)/(nwind-noverlap));	% Number of windows
					% (k = fix(n/nwind) for noverlap=0)

if 0
    disp(sprintf('   x        = (length %g)',length(x)))
    disp(sprintf('   y        = (length %g)',length(y)))
    disp(sprintf('   nfft     = %g',nfft))
    disp(sprintf('   Fs       = %g',Fs))
    disp(sprintf('   window   = (length %g)',length(window)))
    disp(sprintf('   noverlap = %g',noverlap))
    if ~isempty(p)
        disp(sprintf('   p        = %g',p))
    else
        disp('   p        = undefined')
    end
    disp(sprintf('   dflag    = ''%s''',dflag))
    disp('   --------')
    disp(sprintf('   k        = %g',k))
end

index = 1:nwind;
KMU = k*norm(window)^2;	% Normalizing scale factor ==> asymptotically unbiased
% KMU = k*sum(window)^2;% alt. Nrmlzng scale factor ==> peaks are about right

Spec = zeros(nfft,1); Spec2 = zeros(nfft,1);
for i=1:k
    if strcmp(dflag,'none')
        xw = window.*(x(index));
    elseif strcmp(dflag,'linear')
        xw = window.*detrend(x(index));
    else
        xw = window.*detrend(x(index),0);
    end
    index = index + (nwind - noverlap);
    Xx = abs(fft(xw,nfft)).^2;
    Spec = Spec + Xx;
    Spec2 = Spec2 + abs(Xx).^2;
end

% Select first half
if ~any(any(imag(x)~=0)),   % if x is not complex
    if rem(nfft,2),    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
    Spec = Spec(select);
    Spec2 = Spec2(select);
%    Spec = 4*Spec(select);     % double the signal content - essentially
% folding over the negative frequencies onto the positive and adding.
%    Spec2 = 16*Spec2(select);
else
    select = (1:nfft)';
end
freq_vector = (select - 1)*Fs/nfft;

% find confidence interval if needed
if (nargout == 3)|((nargout == 0)&~isempty(p)),
    if isempty(p),
        p = .95;    % default
    end
    % Confidence interval from Kay, p. 76, eqn 4.16:
    % (first column is lower edge of conf int., 2nd col is upper edge)
    confid = Spec*chi2conf(p,k)/KMU;

    if noverlap > 0
        disp('Warning: confidence intervals inaccurate for NOVERLAP > 0.')
    end
end

Spec = Spec*(1/KMU);   % normalize

% set up output parameters
if (nargout == 3),
   Pxx = Spec;
   Pxxc = confid;
   f = freq_vector;
elseif (nargout == 2),
   Pxx = Spec;
   Pxxc = freq_vector;
elseif (nargout == 1),
   Pxx = Spec;
elseif (nargout == 0),
   if ~isempty(p),
       P = [Spec confid];
   else
       P = Spec;
   end
   newplot;
   plot(freq_vector,10*log10(abs(P))), grid on
   xlabel('Frequency'), ylabel('Power Spectrum Magnitude (dB)');
end

function [msg,x,y,nfft,Fs,window,noverlap,p,dflag] = ...
                               psdchk(P1,P2,P3,P4,P5,P6,P7,P8)
%PSDCHK Helper function for PSD, CSD, COHERE, TFE and SPECGRAM
%	assumes that there are between 1 and 8 arguments on input

%	Author(s): T. Krauss, 10-28-93
%	Copyright (c) 1984-94 by The MathWorks, Inc.
%       $Revision: 1.2 $  $Date: 1994/01/25 17:59:39 $

%if 0,   % forces the interpreter to pass parameters by value, not reference
%    P1=P1; P2=P2; P3=P3; P4=P4; P5=P5; P6=P6; P7=P7; P8=P8;
%end

msg = [];
if (nargin>1),
    % 0ne signal or 2 present?
    if (all(size(P1)==size(P2))&~isstr(P2)) | (max(size(P2))>1 & ~isstr(P2))
        onesig = 0;      % two signals, x and y, present
        x = P1; y = P2;
        for k = 3:nargin,    % shift parameters one left
            eval(['P' num2str(k-1) ' = P' num2str(k) ';']);
        end
    else 
        onesig = 1;      % only one signal, x, present
        x = P1; y = []; 
    end
else  % nargin == 1
    onesig = 1;      % only one signal, x, present
    x = P1; y = []; 
end

% now x and y are defined; let's get the rest

if nargin == (1 + (1-onesig)),
    nfft = min(length(x),256);
    window = hanning(nfft);
    noverlap = 0;
    Fs = 2;
    p = [];
    dflag = 'linear';
elseif nargin == (2 + (1-onesig)),
    if isempty(P2),   dflag = 'linear'; nfft = min(length(x),256); 
    elseif isstr(P2), dflag = P2;       nfft = min(length(x),256); 
    else              dflag = 'linear'; nfft = P2;   end
    Fs = 2;
    window = hanning(nfft);
    noverlap = 0;
    p = [];
elseif nargin == (3 + (1-onesig)),
    if isempty(P2), nfft = min(length(x),256); else nfft=P2;     end
    if isempty(P3),   dflag = 'linear'; Fs = 2;
    elseif isstr(P3), dflag = P3;       Fs = 2;
    else              dflag = 'linear'; Fs = P3; end
    window = hanning(nfft);
    noverlap = 0;
    p = [];
elseif nargin == (4 + (1-onesig)),
    if isempty(P2), nfft = min(length(x),256); else nfft=P2;     end
    if isempty(P3), Fs = 2;     else    Fs = P3; end
    if isstr(P4)
        dflag = P4;
        window = hanning(nfft);
    else
        dflag = 'linear';
        window = P4;
        if length(window) == 1, window = hanning(window); end
        if isempty(window), window = hanning(nfft); end
    end
    noverlap = 0;
    p = [];
elseif nargin == (5 + (1-onesig)),
    if isempty(P2), nfft = min(length(x),256); else nfft=P2;     end
    if isempty(P3), Fs = 2;     else    Fs = P3; end
    window = P4;
    if length(window) == 1, window = hanning(window); end
    if isempty(window), window = hanning(nfft); end
    if isstr(P5)
        dflag = P5;
        noverlap = 0;
    else
        dflag = 'linear';
        if isempty(P5), noverlap = 0; else noverlap = P5; end
    end
    p = [];
elseif nargin == (6 + (1-onesig)),
    if isempty(P2), nfft = min(length(x),256); else nfft=P2;     end
    if isempty(P3), Fs = 2;     else    Fs = P3; end
    window = P4;
    if length(window) == 1, window = hanning(window); end
    if isempty(window), window = hanning(nfft); end
    if isempty(P5), noverlap = 0; else noverlap = P5; end
    if isstr(P6)
        dflag = P6;
        p = [];
    else
        dflag = 'linear';
        if isempty(P6), p = .95;    else    p = P6; end
    end
elseif nargin == (7 + (1-onesig)),
    if isempty(P2), nfft = min(length(x),256); else nfft=P2;     end
    if isempty(P3), Fs = 2;     else    Fs = P3; end
    window = P4;
    if length(window) == 1, window = hanning(window); end
    if isempty(window), window = hanning(nfft); end
    if isempty(P5), noverlap = 0; else noverlap = P5; end
    if isempty(P6), p = .95;    else    p = P6; end
    if isstr(P7)
        dflag = P7;
    else
        msg = 'DFLAG parameter must be a string.'; return
    end
end

% NOW do error checking
if (nfft<length(window)), 
    msg = 'Requires window''s length to be no greater than the FFT length.';
end
if (noverlap >= length(window)),
    msg = 'Requires NOVERLAP to be strictly less than the window length.';
end
if (nfft ~= abs(round(nfft)))|(noverlap ~= abs(round(noverlap))),
    msg = 'Requires positive integer values for NFFT and NOVERLAP.';
end
if ~isempty(p),
    if (prod(size(p))>1)|(p(1,1)>1)|(p(1,1)<0),
        msg = 'Requires confidence parameter to be a scalar between 0 and 1.';
    end
end
if min(size(x))~=1,
    msg = 'Requires vector (either row or column) input.';
end
if (min(size(y))~=1)&(~isempty(y)),
    msg = 'Requires vector (either row or column) input.';
end
%if (length(x)~=length(y))&(~isempty(y)),
%    msg = 'Requires X and Y be the same length.';
%end

dflag = lower(dflag);
if strcmp(dflag,'l')|strcmp(dflag,'li')|strcmp(dflag,'lin')| ...
   strcmp(dflag,'line')|strcmp(dflag,'linea')|strcmp(dflag,'linear')
      dflag = 'linear';
elseif strcmp(dflag,'m')|strcmp(dflag,'me')|strcmp(dflag,'mea')| ...
   strcmp(dflag,'mean')
      dflag = 'mean';
elseif strcmp(dflag,'n')|strcmp(dflag,'no')|strcmp(dflag,'non')| ...
   strcmp(dflag,'none')
      dflag = 'none';
else
    msg = 'DFLAG must be ''linear'', ''mean'', or ''none''.';
end

function w = hanning(n_est,sflag)
%HANNING Hanning window.
%   W = HANNING(N) returns the N-point symmetric Hanning window 
%       in a column vector. 
%   W = HANNING(N,SFLAG) generates the N-point Hanning window 
%       using SFLAG window sampling. SFLAG may be either 'symmetric' 
%       or 'periodic'. By default, 'symmetric' window sampling is used. 

%   Copyright (c) 1988-98 by The MathWorks, Inc.
%       $Revision: 1.10 $  $Date: 1997/11/26 20:13:47 $

error(nargchk(1,2,nargin));

[n,w,trivalwin] = check_order(n_est);
if trivalwin, return, end;

% Set sflag to default if it's not already set:
if nargin == 1,
   sflag = 'symmetric';
end

switch lower(sflag),
case 'periodic'
   w = sym_hanning(n+1);
   w = w(1:end-1);
case 'symmetric'
   w = sym_hanning(n);
otherwise
	error('Sampling must be either ''symmetric'' or ''periodic''');
end

function w = sym_hanning(n)
w = .5*(1 - cos(2*pi*(1:n)'/(n+1)));

function [n_out, w, trivalwin] = check_order(n_in)
%CHECK_ORDER Checks the order passed to the window functions.
% [N,W,TRIVALWIN] = CHECK_ORDER(N_ESTIMATE) will round N_ESTIMATE to the
% nearest integer if it is not alreay an integer. In special cases (N is [],
% 0, or 1), TRIVALWIN will be set to flag that W has been modified.

%   Copyright (c) 1988-98 by The MathWorks, Inc.
%   $Revision: 1.5 $  $Date: 1997/11/26 20:14:55 $

w = [];
trivalwin = 0;

% Special case of negative orders:
if n_in < 0,
   error('Order cannot be less than zero.');
end

% Check if order is already an integer or empty
% If not, round to nearest integer.
if isempty(n_in) | n_in == floor(n_in),
   n_out = n_in;
else
   n_out = round(n_in);
   warning('Rounding order to nearest integer.');
end

% Special cases:
if isempty(n_out) | n_out == 0,
   w = zeros(0,1);               % Empty matrix: 0-by-1
   trivalwin = 1; 
elseif n_out == 1,
   w = 1;
   trivalwin = 1;   
end
